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Abstract 

Fiber strength is the key trait that determines fiber quality in cotton, and it is closely related to secondary cell wall synthesis. 
To understand the mechanism underlying fiber strength, we compared fiber transcriptomes from different G. barbadense 
chromosome introgression lines (CSILs) that had higher fiber strengths than their recipient, G. hirsutum acc. TM-1. A total of 
18,288 differentially expressed genes (DEGs) were detected between CSIL-35431 and CSIL-31010, two CSILs with stronger 
fiber and TM-1 during secondary cell wall synthesis. Functional classification and enrichment analysis revealed that these 
DEGs were enriched for secondary cell wall biogenesis, glucuronoxylan biosynthesis, cellulose biosynthesis, sugar-mediated 
signaling pathways, and fatty acid biosynthesis. Pathway analysis showed that these DEGs participated in starch and sucrose 
metabolism (328 genes), glycolysis/gluconeogenesis (122 genes), phenylpropanoid biosynthesis (101 genes), and oxidative 
phosphorylation (87 genes), etc. Moreover, the expression of MYB- and NAC-type transcription factor genes were also 
dramatically different between the CSILs and TM-1. Being different to those of CSIL-31 134, CSIL-35431 and CSIL-31010, there 
were many genes for fatty acid degradation and biosynthesis, and also for carbohydrate metabolism that were down- 
regulated in CSIL-35368. Metabolic pathway analysis in the CSILs showed that different pathways were changed, and some 
changes at the same developmental stage in some pathways. Our results extended our understanding that carbonhydrate 
metabolic pathway and secondary cell wall biosynthesis can affect the fiber strength and suggested more genes and/or 
pathways be related to complex fiber strength formation process. 
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Introduction 

The cotton fiber is a terminally differentiated single cell derived 
from the epidermal cell of the developing ovule. After initiation, 
the fiber cell undergoes 1000- to 3000-fold elongation during its 
development. The development of cotton fibers involves four 
partially overlapping stages: initiation (—3 to +3 days post- 
anthesis; DPA), elongation and primary cell wall formation (3—23 
DPA), secondary cell wall formation (16-40 DPA) and maturation 
(40-50 DPA) [1-6]. The most rapid period of fiber cell elongation 
begins around 10-16 DPA and continues to ~20 DPA. Primary 
and secondary cell wall synthesis overlaps during the period of 1 6- 
25 DPA. During the secondary cell wall formation stage, the speed 
of cell elongation slows down and even stops. 

Fiber strength is an important indicator of cotton fiber quality, 
and depends on formation of the secondary cell wall. Cellulose 
synthesis plays a predominant role in fiber cells, and cellulose 
accounts for >95% of the dry weight of the mature cotton fiber 



[3,7] . Genome and EST sequencing have revealed that there are 
at least ten different CesA genes for cellulose synthase in 
Arabidopsis; CesA-like genes have also been reported in rice and 
barley [8-10]. In cotton [Gossypium raimondii), at least 15 cellulose 
synthase (CESA) sequences are required for cellulose synthesis 
[1 1]. A recent investigation in Arabidopsis thaliana using microarrays 
led to the identification of genes that are highly co-expressed with 
cellulose synthase genes and two mutants, irx8 and irxl3, that 
have irregular xylem phenotypes, were also identified [12]. 
Sucrose synthase (Susy) is the enzyme that catalyzes the hydrolysis 
of sucrose to UDP-glucose that is then used as a substrate for 
cellulose synthesis. In cotton, the expression of Susy is higher at 
16-32 DPA, and this enzyme plays a major role in partitioning 
carbon toward cellulose synthesis in the fiber [13]. SusC is another 
new sucrose synthase gene with a high level of expression during 
secondary cell wall synthesis [14], Peroxide, mainly as H 2 0 2 , 
promotes cellulose synthesis as a signal of secondary cell wall 
synthesis [15,16]. 
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At present, many ovule- and fiber-specific cDNA libraries have 
been constructed and sequenced, and more than 268,000 
expressed sequence tags (ESTs) from Gossypium are deposited in 
the NCBI database (http://www.ncbi.nlm.nih.gov). For genetic 
characterization of rapid cell elongation in cotton fibers, approx- 
imately 14,000 unique genes were assembled from 46,603 
expressed sequence tags (ESTs) from developmentally-staged fiber 
cDNAs of a cultivated diploid species (G. arboreum L.). Eighty-one 
genes that were significantly up-regulated during secondary cell 
wall synthesis were found to be involved in cell wall biogenesis and 
energy/carbohydrate metabolism, which is consistent with the 
stage of cellulose synthesis during secondary cell wall modification 
in developing fibers [17]. Transcriptome profiling of the cotton 
fiber early in development by high-throughput tag-sequencing 
(Tag-seq) analysis using the Solexa Genome Analyzer reveals 
significant differential expression of genes in a fuzzless/lindess 
mutant [18]. High-throughput, genome- wide transcriptomic 
analysis of cotton under drought stress revealed a significant 
down-regulation of genes and pathways involved in fiber 
elongation, and an up-regulation of defense response genes [19]. 
More research have been processed in fiber initiation and 
elongation stage [20-24]. Saturated very-long-chain fatty acids 
(VLCFAs; G20:0-G30:0) exogenously applied in ovule culture 
medium significantly promoted fiber cell elongation in cotton (G. 
hirsutum L.) by activating ethylene biosynthesis [25,26]. Previous 
investigations into cotton fiber development mainly focused on the 
elongation stage, and the number of genes reported from the later 
stages is quite small. Most of the genes up-regulated during 
secondary cell wall synthesis were related to cellulose synthesis, cell 
wall biosynthesis, and carbohydrate metabolism [17,22,27]. 

Chromosome segment introgression lines (CSILs) consist of a 
battery of near-isogenic lines that have been developed to cover 
the entire genomes of some crops, including tomato, rice, wheat, 
and cotton [28-31]. With the exception of a single, homozygous 
chromosome segment transferred from a donor parent, the 
remaining genome of each CSIL is the same as the recipient 
parent [31]. We used G. barbadense CSILs in the background of the 
standard genetic line of G. hirsutum, cv. TM-1, in order to 
understand the molecular mechanism behind superior quality 
fiber formation. Multi-point tests showed that three CSILs 
produced stronger fibers when compared to the recipient parent 
TM-1, but one CSIL produced weaker fibers. Using Solexa 
Genome sequencing, we analyzed transcriptome profiles from the 
CSILs and TM-1. We found that many genes were either up- or 
down-regulated at the stage of secondary cell wall synthesis, and 
that many metabolic pathways were altered in the CSILs. 

Materials and Methods 

Plant materials 

G. hirsutum cv. TM-1, the genetic standard for Upland cotton, 
was obtained from the Southern Plains Agricultural Research 
Center, USDA-ARS, College Station/Texas, USA [32]. G. 
barbadense cv. Hai7124, an extra-long staple cotton that is widely 
grown in China, is descended from a selected individual in a study 
of inheritance of resistance to Verticillium dahlia [33,34]. In this 
study, we identified three CSILs with stronger fiber or high fiber 
strength that carried different G. barbadense chromosome segment(s) 
in the recurrent parent TM- 1 . The detailed method of developing 
CSILs has been described previously [31]. We selected three 
CSILs, CSIL-35431, CSIL-31 134, and CSIL-31010, in which the 
average fiber strength were 35.1, 34.73 and 34.28 cN/tex, 
respectively, significandy higher than TM-1, and also CSIL- 
35368 which had poorer fiber strength than TM-1 (28. 71 cN/tex) 



(Table SI). The introgressed G. barbadense chromosomal segments 
were different in the four lines [35]. Fiber samples were collected 
at 15, 20, and 25 DPA, frozen in liquid nitrogen, and stored at — 
70°C. 

RNA isolation and evaluation 

Total RNA was extracted from frozen tissue using an improved 
CTAB extraction protocol [36] . RNAs were evaluated for quality 
using RNA Pico Chips on an Agilent 2100 Bioanalyzer (Agilent 
Technologies, Santa Clara, CA, USA). All RNA samples were 
quantified and qualified with an RNA Integrity Number (RIN) > 
8, and 28S/18S rRNA band intensity (2:1). 

Library construction and sequencing 

Digital gene expression libraries were constructed using the 
Illumina Gene Expression Sample Preparation Kit according to 
the manufacturer's instructions. We constructed and sequenced 14 
libraries derived from immature fibers at 15, 20, and 25 DPA 
using the Solexa Genome Sequencing Analyzer system provided 
by BGI (Beijing Genomics Institute at Shenzhen, China), which 
gave 21 bp tags. The process was described in detail previously 
[18]. 

Data processing, statistical evaluation, and selection of 
differentially expressed genes (DEGs) 

Raw data reads were filtered by the Illumina pipeline to 
produce clean data. All low-quality data, such as short tags (< 
21 nt) and singletons, were removed. A database of 2 1 -base-long 
sequences was produced beginning with CATG using 37,505 
reference genes from the diploid species G. raimondii (http:/ /www. 
phytozome.net). The remaining high quality sequences were then 
mapped to this database; only a single mismatch was allowed, and 
more than one match was excluded. Gene expression levels were 
the summation of tags aligned to the different positions of the same 
gene. Expression levels are expressed as TPM, transcripts per 
million. To identify DEGs during fiber elongation, we compared 
pairs of DEG profiles from different libraries. Three fiber 
development periods for the four CSILs were compared with 
the same period for TM- 1 , and 1 1 comparisons were obtained. P- 
and Qjvalues were also calculated for every comparison [37]. 
DEGs were defined as FDR<0.001 with an absolute value of 
| log 2 Ratio | & 1 to judge the significance of differences in transcript 
abundance. 

Digital tag profiling analysis 

DEG clustering in CSILs at different developmental stages were 
performed with Cluster3.0 (http://bonsai.hgc.jp/~mdehoon/ 
software/cluster/software.htm). We also performed clustering 
with the 'Self-organizing tree algorithm' (SOTA, Multiple Array 
Viewer software, MeV 4.9.0) [38]. 

GO enrichment and KEGG (Kyoto Encyclopedia of Genes and 
Genomes) pathway analysis was done using BLAST2GO (http:// 
www.blast2go.com/b2ghome). Mapman was also used to analyze 
metabolic pathway base on KEGG database [39]. 

Quantitative RT-PCR 

Quantitative RT-PCR assays were performed on a 7500 Real- 
Time PCR system (Applied Biosystems, San Francisco, CA, USA). 
Reactions were performed in a final volume of 15 uL and 
contained 2 uL of diluted cDNA, 7.5 uL of 2 x SYBR mix (Roche, 
Basel, Switzerland), and 200 nM of the forward and reverse 
primers. Primer lengths were designed to be between 18 and 24 nt 
using Beacon Designer 7, and PCR amplicon lengths were 
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designed to be between 100 bp and 150 bp (Table S2). The 
thermal cycling conditions were 40 cycles of 95°C for 15 s, 60°C 
for 30 s, and 72°C for 30 s. All reactions were run in triplicate, 
and the cotton histoneS gene (ACC NO. AF024716) was used as an 
internal control for normalization of expression levels (F: 5'- 
GGTGGTGTGAAGAAGCCTCAT-3 ', and R: 5'-AATTT- 
CACGAACAAGCCTCTGGAA-3'). The relative gene expres- 
sion levels were presented as 2 ACT . 

Results 

Statistical analysis of transcriptome data 

The total number of sequence tags per library ranged from 7.0 
to 8.5 million, and the number of distinct sequence tags was 
between 1.8 and 2.2 million. Approximately 50% of the clean tags 
were mapped to reference genes, and 60 % of the reference genes 
were mapped with unambiguous tag (Table 1 and Table S3). 

To see whether the fiber transcriptomes at different develop- 
mental stages were different, the 23,237 genes which were 
expressed in at least three libraries at one stage (15 DPA, 20 
DPA, or 25 DPA) were classified into six groups using the Multiple 
Array Viewer using TPM value (Figure 1A). Genes in Group 3 
had higher expression levels at 15 DPA and 20 DPA than at the 
later stage (25 DPA). Genes in Group 4 had higher expression 
levels at 15 DPA than at either 20 DPA or 25 DPA. Genes in 
Group 5 showed the opposite expression pattern, with higher 
expression levels at 20 DPA and 25 DPA compared to 15 DPA. 
The other groups also showed distinct expression patterns 
(Figure 1A). 

Classification by gene function revealed that Group 3 is 
enriched in genes involved in protein catabolism, cell division, 
and cellulose biosynthesis, Group 4 is enriched in genes for cell 
morphogenesis, fatty acid biosynthesis, lipid transport, and wax 
biosynthesis, and Group 5 has more genes involved in glucose 
catabolism, response to chitin, and sucrose metabolism (Figure IB). 
The unbalanced pattern of the expressed-gene functional distri- 
bution could possibly reflect some physiological events involved in 
secondary cell wall biosynthesis. 

Cluster analysis of differentially expressed genes (DEGs) 
between and/or among CSILs 

We specifically looked for DEGs in secondary cell wall fibers 
from 15 to 25 DPA, because previous studies have reported that 
the different sets of transcripts responsible for fiber secondary cell 
wall formation may be enriched at these stages of development 
[17,22,27]. Three fiber development periods for the four CSILs 
were compared with TM-1 at the same period. DEGs were 
defined as FDR<0.001 with an absolute value of |log2Ratio| Si. 
Analysis of the data indicated that many genes showed differential 
expression in the 1 1 comparison groups. The number of DEGs 
were about 6,000-8,000 in CSILs from 15 DPA to 25 DPA 
(Figure 2A). But the number of DEGs in CSIL-31010 at 20 DPA, 
CSIL-31010 at 25 DPA, and CSIL-31 134 at 15 DPA, were 4,600, 
10,106 and 2,060, respectively. We also found that the DEGs that 
were up-regulated or down-regulated were different in CSILs. 
There were -1,500-3,500 DEGs in common from 15 DPA to 25 
DPA between CSIL-35431, CSIL-31010, and CSIL-35368 
(Figure 2B). 

To understand the mechanisms behind the changes in fiber 
strength observed in the CSILs, we also analyzed the common 
DEGs among CSIL-35431, CSIL-31010 and CSIL-31 134 (Table 
S4). A total of 727 and 1796 common DEGs were selected at 15 
and 20 DPA in three stronger fiber CSILs, respectively (Figure 3). 
More functional enrichment were shown at 15 DPA, including 
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Figure 1 . Statistical analysis of transcriptome data. (A) SOTA clustering of the different genes using Log2(TPM). T, TM-1 ; A, CSIL-35431 ; B, CSIL- 
31010; C, CSIL-31 134; D, CSIL-35368. 15, 15 DPA; 20, 20 DPA; 25, 25 DPA. (B) Distribution of functions of genes in different clusters. Yellow square 
indicated group 3, green square indicated group 4 and blue square indicated group 5. X-axis indicated different enriched process and Y-axis indicated 
number of hit-found genes in these processes. 
doi:1 0.1 371 /journal.pone.0094642.g001 



major CHO metabolism (carbohydrate), cell wall biosynthesis, 
amino acid metabolism and secondary metabolism (Figure 3E). 
Among these genes, 321 and 998 common upregulated DEGs 
between the same CSILs at 15 and 20 DPA were indentified, 
respectively (Figure 3). These common DEGs or processes maybe 
directly related to the fiber strength. However, these DEGs maybe 
function as downstream genes altered by the introgressed segments 
since these CSILs were inserted different G. barbadmse segments in 
recipient TM- 1 . 

To visualize the expression patterns of DEGs, we performed 
cluster analysis of 18,288 genes that were differentally expressed 
between CSIL-35431 and CSIL-31010 (Figure 4). These DEGs 
could be grouped into six clusters, designated G1-G6, based on 
their expression patterns. From 15 DPA to 20 DPA, the stages of 
fast fiber elongation and secondary cell wall deposition overlap, 
with the latter reaching a peak at around 20-25 DPA. We focused 
on clusters Gl, G4, and G6 to conduct data analysis in order to 
identify genes that were either up-regulated or down-regulated 
during the secondary cell wall synthesis stage. Compared to the 
TM-1 control, 3,658 genes in cluster Gl were highly expressed at 
15 and 20 DPA, 4,487 genes in G4 were highly expressed at 15 
DPA, 20 DPA, and 25 DPA, 3,033 genes in G6 were highly 
expressed only at 25 DPA, and the other three groups showed 
various different expression patterns. Clustering results for 19,742 



DEGs from the four CSILs showed five groups, indicating that the 
gene expression pattern in CSIL-31 134 was distinct from the 
others at 15 DPA and 20 DPA, and that CSIL-35368 was similar 
to CSIL35431 and CSIL-31010 (Figure SI). 

Functional annotation by GO enrichment and KEGG 
analysis 

To understand the mechanisms behind the changes in fiber 
strength observed in the CSILs, we analyzed DEG enrichment in 
the major functional GO categories of biological process, 
molecular function, and cellular component between CSIL- 
35431 and CSIL-31010. Based on the clustering results shown 
in Figure 4, Gl was enriched in genes for secondary cell wall 
biogenesis, glucuronoxylan biosynthesis, microtubule-based move- 
ment, and cellulose biosynthesis, G4 was enriched in genes for 
protein phosphorylation, response to chitin, and sugar-mediated 
signaling pathways, and G6 was enriched in fatty acid biosynthesis 
genes (Table 2). These data suggest that in the developmental 
stage of secondary cell wall deposition, DEGs were enriched for 
carbohydrate synthesis and cell wall formation. 

We applied the same GO analysis to the common DEGs at 15 
DPA and 20 DPA in CSIL-35431 and CSIL-31010, respectively. 
These DEGs were enriched in genes for similar functional 
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Figure 2. Statistical of DEGs between CSILs and TM-1 at 15, 20 and 25 DPA. (A) Up-regulated and down-regulated genes in different 
comparison. Red bar, up-regulated genes compared to TM-1; green bar, stand for down-regulated genes compared to TM-1, Blue square, total DEGs. 
CSILs included CSIL-35431, CSIL-3010, CSIL-31134, CSIL-35368 and TM-1. 15, 15DPA; 20, 20DPA; 25, 25DPA. (B) Common and special DEGs at 15 DPA, 
20 DPA and 25 DPA. 
doi:1 0.1 371 /journal.pone.0094642.g002 



categories, such as cellular metabolic processes and carbohydrate 
metabolism, etc. We also found genes for some processes that were 
enriched only in CSIL35431 or CSIL-31010 (Figure S2). 

Further GO analysis for CSIL-35368 and CSIL-31 134 indicat- 
ed that the DEGs in CSIL-35368 at 15 and 20 DPA were enriched 
in genes for lignin biosynthesis, secondary cell wall biogenesis, and 
response to chitin, which was similar to the enrichment found in 
CSIL-35431 and CSIL-31010. But at 15 and 20 DPA in the 
stronger fiber line CSIL-31134, GO enrichments were different 
from the other three lines, mainly in genes for ATP synthesis, 



proton transport, copper ion export, and oxidoreductase activity, 
but not in cell wall biosynthesis (Table S5). 

Based on the results of GO analysis, we know that the secondary 
cell wall related biological process were impacted in the CSILs, but 
it is still not very clear how secondary cell wall biosynthesis was 
affected in the CSILs. Therefore, we performed pathway analysis 
on 18,288 DEGs in CSIL-35431 and CSIL-31010. The most 
highly enriched pathways found are listed in Table 3. KEGG 
analysis showed that the genes were enriched in pathways for 
starch and sucrose metabolism (328 genes), glycolysis/gluconeo- 
genesis (122 genes), phenylpropanoid biosynthesis (101 genes), and 
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Figure 3. Analysis of common and common upregulated DEGs among three stronger fiber CSILs. (A, B, C, D) Common and common 
upregulated DEGs among three stronger fiber CSILs at 1 5 and 20 DPA. Common_up, common regulated DEGs. (B) Functional enrichment analysis of 
these DEGs using mapman software (Summary statistic type, wlcoxon). Colors from blue to red indicated that functions were enriched more 
significantly with smaller p-values. 
doi:1 0.1 371 /journal.pone.0094642.g003 



oxidative phosphorylation (87 genes) (Table 3 and Figure S3). The 
regulation of some enzymes that catalyze sucrose, starch, and 
cellulose biosynthesis may have a direct or indirect impact on fiber 
quality. This could be especially true for sucrose and pectin 
metabolism, and many genes in these pathways were up-regulated. 



We also found that genes involved in phenylpropanoid and 
flavonoid biosynthetic processes were enriched in the CSILs. 

Based on the cluster analysis of the weaker fiber line CSIL- 
35368, we hypothesized that changes in other biochemical 
pathways led to reduced fiber strength (Figure SI). Considering 
only those that were down-regulated in CSIL-35368, we found 
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Figure 4. Heat map analysis of the expression of DEGs between CSILs and TM-1. A, B and T indicated CSIL-35431, CSIL-31010 and TM-1, 
respectively. 1 5, 1 5DPA; 20, 20DPA; 25, 25DPA. Red color indicated up-regulated genes and green color indicated down-regulated genes. N = number 
of DEGs in different group. 
doi:1 0.1 371 /journal.pone.0094642.g004 
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Table 2. Enrichment analysis of gene ontologies from 15 to 25 DPA. 



Cluster 


GO-ID 


GO Ontology (Biological process) 


G1 


GO:0010417 


glucuronoxylan biosynthetic process 


15DPA up-regulated 


GO:0009834 


secondary cell wall biogenesis 


20DPA up-regulated 


GO:0007018 


microtubule-based movement 


25DPA down-regulated 


GO:0030244 


cellulose biosynthetic process 




GO:0009753 


response to jasmonic acid stimulus 


G2 


GO:0015031 


protein transport 


15DPA down-regulated 


GO:0015991 


ATP hydrolysis 


20DPA down-regulated 


GO:0032544 


plastid translation 


25DPA down-regulated 


GO:0016075 


rRNA catabolic process 




GO:000651 1 


ubiquitin catabolic process 


G3 


GO:0009734 


auxin mediated signaling pathway 


15DPA down-regulated 


GO:0030259 


lipid glycosylation 


20DPA up-regulated 


GO:0018106 


peptidyl-histidine phosphorylation 


25DPA down-regulated 


GO:0009722 


detection of cytokinin stimulus 


G4 


GO:0006468 


protein phosphorylation 


15DPA up-regulated 


GO:0010200 


response to chitin 


20DPA up-regulated 


GO:0006096 


glycolysis 


25DPA up-regulated 


GO:0010182 


sugar mediated signaling pathway 




GO:0009966 


regulation of signal transduction 




LjU.UU 1 Doo^f 


protein folding 


15DPA up-regulated 


GO:0006886 


intracellular protein transport 


20DPA down-regulated 


GO:0006122 


mitochondrial electron transport 


25DPA up-regulated 


GO:0015914 


phospholipid transport 


G6 


GO:0007267 


cell-cell signaling 


15DPA down-regulated 


GO:0010025 


wax biosynthetic process 


20DPA down-regulated 


GO:0006633 


fatty acid biosynthetic process 


25DPA up-regulated 


GO:0006723 


hydrocarbon biosynthetic process 



G1-G6 according to Figure 3. 

doi:1 0.1 371 /journal.pone.0094642.t002 



genes that participated in fatty acid degradation and biosynthesis, 
and also in carbohydrate metabolic pathways (Figure 2B and 
Figure S4). 

Eight genes previously reported in the carbohydrate pathway 
were selected for quantitative RT-PCR. The expression patterns 
of these genes were consistent with the DEG data in TM-1 
(Figure 5) and in the CSILs as well (Figure 6B and Figure 7B). 

Carbohydrate metabolism in the secondary cell wall 
synthesis stage 

Following the start of secondary cell wall formation, protein and 
carbohydrate metabolism genes involved in cell wall biosynthesis 
will be up-regulated [27]. We selected 72 DEGs associated with 
carbohydrate metabolism to investigate the mechanism of fiber 
development. These genes were related to pectin, sucrose, 
galactan, glucan, xyloglucan, and cellulose biosynthesis. We were 
interested in genes that are up-regulated in fiber cells at 15 DPA 
and 20 DPA, at the start of secondary cell wall formation. A heat 
map showing the different expression levels for these genes 
including cellulose synthase, sucrose synthase, pectin lyase, and 
other polysaccharides degradation in CSIL-35431 and CSIL- 
3 1 0 1 0 is shown in Figure 6A. We found that the cellulose synthase 
genes were up-regulated in the CSILs at 15 DPA-25 DPA. It has 



been reported that cellulose biosynthesis predominates, and that 
many other metabolic pathways are down-regulated during 
secondary cell wall synthesis [27]. Moreover, we confirmed the 
expression patterns of cellulose synthase genes, annotated with the 
Arabidopsis genes AtCESA4, AtCESA7 and AtCESAS, using quanti- 
tative RT-PCR (Figure 4B). Proteins encoded by AtCESA4, 7, and 
8 are specifically required to form a functional cellulose synthase 
complex (CSC) that is essential for secondary cell wall formation 
[40-42]. 

Transcription factors associated with secondary cell wall 
synthesis 

Recent molecular and genetic studies have identified transcrip- 
tion factors that are involved in regulating secondary cell wall 
synthesis in Arabidopsis [43-45]. In our study, 97 MYB-type and 68 
NAC-type transcription factors showed changes in expression 
between the CSILs and TM-1 (Table S6, Table S7). It was 
interesting that some NACs and MYBs were up-regulated in 
CSIL-35431 and CSIL-31010 during the secondary ceU wall 
synthesis stage, especially at 15 DPA and 20 DPA. Defined as 
|log 2 Ratio| >2, 59 MYB and 47 NAC transcription factors were 
selected for heat-map analysis (Figure 7A). Among these 
transcription factors, genes homologous with ATMTB2, AT- 
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Figure 5. Quantitative RT-PCR validation of tag-mapped genes in TM-1. These genes have been reported before, including 3 CesA genes 
(A,B,C) (homologous with AtCESA4, AtCESA7, AtCESA8, respectively), xyloglucan endotransglucosylase (D), beta -galactosidase (E), glycosyl hydrolase 
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MTB43, ATMYB73, ATNAC52, and ATMAC61 were expressed at 
higher levels in the CSILs. We confirmed that three transcription 
factors were up-regulated in CSILs from 15 DPA to 25 DPA 
(Figure 7B). In the MYB family, it has been reported that the 
expression of genes for MYB85, MYB52, MYB54, MYB69, 
MYB42, and MYB43 are developmentally associated with cells 
undergoing secondary wall thickening [45] . 

Different metabolic pathways associated with altered 
fiber strength 

In order to investigate the mechanisms underlying changes in 
fiber strength, we analyzed several metabolic pathways including 
cell wall, lipids, minor CHO (carbohydrate) metabolism, and two 
secondary metabolite pathways. It is interesting that DEGs 



involved in cell wall proteins, cell wall pectin esterase, cell wall 
modification, cell wall cellulose synthesis, cell wall degradation/ 
pectate lyases, lipid metabolism/FA synthesis, and lipid degrada- 
tion showed distinct expression patterns or differential up/down- 
regulation at 20 DPA (Figure 8A). We found that up-regulated 
DEGs were similar to down-regulated DEGs both in CSIL-35431 
and CSIL-35368. However, most of DEGs in CSIL-31010 were 
up-regulated at 20 DPA, while the opposite was true for DEGs in 
CSIL-31134, especially those genes involved in cell wall modifi- 
cation. In CSIL-31134, we also found a few genes in these 
metabolic pathways that were changed at 15 DPA except in cell 
wall modification, and in CSIL-31010, we found DEGs enriched 
in these metabolic pathways at 25 DPA (Figure S5). From the 
secondary metabolism results, we identified a few DEGs involved 
in flavonoid biosynthesis in CSIL-35431 and CSIL-31010 at 15 
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DPA. In contrast, more genes were up-regulated or down- 
regulated in CSIL-35368 at 15 DPA. It was obvious that DEGs 
from the phenylpropanoid pathways at 25 DPA were different 
from one another, and the expression pattern of DEGs in CSIL- 
31010 changed dramatically. Moreover, there were few genes that 
were up-regulated or down-regulated in CSIL-35368 at 25 DPA 
(Figure 7B). We assume that metabolic pathways in the CSILs at 
different developmental stages were changed in various ways as a 
result of the introgressed chromosmal segments from G. barbademe. 

Discussion 

G. hirsutum produces a high yield of cotton with moderate fiber 
strength. G. barbadense is characterized by a low yield, but with 
increased fiber fineness and strength. As a breeding target, we tried 
to combine the high yield of G. hirsutum with the superior fiber 



qualities of G. barbadense, and we also wanted to elucidate the 
molecular mechanism behind the formation of superior quality 
fibers. Fiber strength is an important indicator of the cotton fiber 
quality, and depends on the formation of the secondary cell wall. 
Genome-wide transcriptome profiling is effective at revealing 
significant genes and pathways involved in secondary cell wall 
formation. Transcriptome analysis showed that gene expression 
patterns and functional distribution were different during second- 
ary cell wall biosynthesis. 

Carbohydrate metabolism plays an important role in 
secondary cell wall synthesis 

It is well known that the mature cotton fiber is composed of 
nearly pure cellulose, and that such a high level of cellulose 
synthesis requires an abundant supply of UDP-glucose [46,47]. 
This means that a large amount of cellulose is required during the 
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secondary cell wall synthesis stage. Functional classification and 
enrichment analysis showed that following the initiation of 
secondary cell wall synthesis, DEGs were enriched for secondary 
cell wall biogenesis, glucuronoxylan biological processes, and other 
carbohydrate metabolic pathways in the CSILs (Table 2). 
Focusing on carbohydrate metabolic pathways, it is obvious that 
the key intermediate in the multiple pathways is UDP-glucose, a 
substrate for cellulose synthesis. Our results showed that several 
CesA genes are expressed at higher levels during secondary cell 
wall synthesis than they are at earlier stages (Figure 6B). Ten 
AtCESA genes have been reported in Arabidopsis, and AtCESA4, 7, 
and 8 are specifically required to form the cellulose synthase 
complex (CSC) that is essential for secondary cell wall synthesis 
[40-42]. Similarly, three CESA isoforms have been identified 
during secondary cell wall synthesis in rice, maize, and Populus 
[10,48,49]. Also, many genes that participate in the degradation of 
poly- and oligo-saccharides were found to be up-regulated at 15 
and 20 DPA, in order to produce more UDP-glucose for cellulose 
biosynthesis. Similarly, it has also been reported that during the 
secondary cell wall synthesis stage, certain metabolic pathways, 
including hydrolysis of fatty acids and non-cellulose poly- and 
oligo-saccharides, would be up-regulated [27]. Sucrose synthase 
(SuSy) has long been studied as a cytoplasmic enzyme in plant 
cells, where it serves to degrade sucrose and provide carbon for 



respiration and synthesis of cell wall polysaccharides and starch 
[50]. It has also been shown that genes associated with secondary 
cell wall biosynthesis are involved in sugar metabolism [51]. 

Multiple mechanisms affect fiber strength development 

Except for carbohydrate metabolism, recent research has shown 
that transcription factors also affect fiber development during 
secondary cell wall biosynthesis. Several NAC- and MYB-type 
transcription factors were up-regulated in the CSILs compared to 
TM-1 from 15 DPA to 25 DPA, and these included cotton 
homologs of AtMTB2, AtMW43, and AWAC52 etc. (Figure 7A). 
The NAC-mediated transcriptional regulation of secondary wall 
biosynthesis is a conserved mechanism throughout vascular plants 
[44,52]. SMD2, a NAC transcription factor gene, regulates genes 
involved in secondary cell wall development in Arabidopsis fibers 
and increases fiber cell area in Eucalyptus [53]. A MYB75- 
associated protein complex is likely to be involved in modulating 
secondary cell wall biosynthesis in both the Arabidopsis inflores- 
cence and stem [54] . It has also been found that the rice and maize 
MYB transcription factors, OsMYB46 and ZmMYB46, are 
functional orthologs of Arabidopsis MYB46/MYB83 and, when 
overexpressed in Arabidopsis, are able to activate the entire 
secondary wall biosynthetic program [55]. 
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Several metabolic pathways were examined to determine the 
mechanism behind changes in fiber strength; these included cell 
wall, lipids, minor CHO metabolism, and two secondary 
metabolic pathways. Although results of the GO and KEGG 



analyses showed that CSIL-35431, CSIL-31010, and CSIL-35368 
had similar patterns, fiber strength in these three lines were 
different. Our results support the hypothesis that different 
metabolic pathways can affect fiber strength, and the same 
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pathway in the CSILs can be altered differentially at various times 
in development. DEGs in CSIL-31010 were up-regulated at 20 
DPA, while the opposite was found for DEGs in CSIL-31134, 
especially those genes involved in cell wall modification. The 
expression levels of genes involved in flavonoid biosynthesis in the 
weak fiber line CSIL-35368 were changed dramatically at 15 
DPA, but there were few genes changed at 25 DPA; this patter was 
the opposite of that in CSIL-35431 and CSIL-31010, lines with 
high quality fiber. We hypothesize that phenylpropanoid and 
flavonoid metabolism generally affected the fiber strength of 
CSIL-35368. Genes for phenylpropanoid and flavonoid biosyn- 
thesis showed significant enrichment and temporal differences in 
gene expression patterns which are associated with xylem 
formation [56]. It has been reported that expression levels of 
phenylpropanoid genes showed high correlations with specific 
fiber properties, supporting a role in determining fiber strength 
[57]. 

In conclusion, upland cotton CSILs carrying distinct G. 
barbadense chromosomal segments provide valuable material for 
research into fiber development. The G. barbadense chromosome 
segments resulted in different patterns of differentially expressed 
genes, and altered different metabolic pathways, mainly in 
carbohydrate metabolism. In addition, several transcription factor 
genes were found to be specifically up-regulated in the CSILs. 
Metabolic pathways involved in cell wall, lipid, phenylpropanoid, 
and flavonoid biosynthesis play a significant role during secondary 
cell wall formation, and are associated with the development of 
cotton fiber strength. 

Supporting Information 

Figure SI Heat map of the expression of DEGs between 
4 CSILs at 15-25 DPA. 

(TIF) 

Figure S2 Enrichment analysis of common DEGs at 
15DPA and 20DPA in CSIL-35431 and CSIL-31010. 

(TIF) 
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